Bogota car accidents 2014

Bogota car accidents 2014

library(tidyverse)
library(lubridate)

Data import

Data downloaded from https://docs.google.com/spreadsheets/d/1lgR-UpR9v9lmtoeK1ZhlVfs0FBj1xXz_eImhcrSI4Q0/

accidents_raw <- read_csv("data/Car_Accidents_2014_Bogota.csv")

Parsed with column specification:
cols(
  the_geom = col_character(),
  cartodb_id = col_integer(),
  codigo_accidente = col_integer(),
  fecha_ocurrencia = col_date(format = ""),
  hora_ocurrencia = col_character(),
  gravedad = col_character(),
  objectid = col_integer(),
  localidad = col_character(),
  municipio = col_character(),
  x = col_double(),
  y = col_double(),
  timedate = col_datetime(format = ""),
  clase = col_character()
)

accidents_raw

# Odd to have exactly 10000 rows, if the dataset is supposed to be complete

• the_geom : geom
• cartodb_id : reference#
• codigo_accidente : Accident code from police ?
• fecha_ocurrencia eg. 2014-10-29 : Date
• hora_occurencia eg. 14:10:00 : Time
• gravedad : Gravity (severity)
• objectid : reference#
• localidad : Location
• municipio : Municipality
• x : Longitude
• y : Latitude
• timedate <S3: POSIXct> : eg. 2014-10-29 14:10:00
• clase : Class [choque = shock; atropello = ; otro = ; caida de ocupante = ]

# The data is for one year (2014). If we want to filter by month (etc), need a `month` variable first. Choose to use **lubridate** (because we have `timedate`). Could instead `separate()` the `fecha_ocurrencia` (and) and `hora_occurencia` (time) fields into components, but there are 9 rows where `hora_occurencia` says "1899-12-30" (instead of 00:00:00, that `timedate` says), and we'd have to filter/edit those
accidents <- accidents_raw %>%
    filter(year(timedate) == 2014)            # 18 rows from 2013 filtered out
accidents                                     # (9982 rows x 13 cols)

Basic EDA

# Just to get a sense of what's going on
accidents %>% count(the_geom)            # 6924 rows (not unique, some accidents at same x&y)

accidents %>% count(the_geom) %>% count(n)  # 21 rows (1 to 22, without 18)

accidents %>% count(cartodb_id)          # 9982 rows (unique)

accidents %>% count(codigo_accidente)    # 9982 rows (unique). ~5e+05 (?) or 4e+06 (`objectid`)

accidents %>% count(month(timedate))     # 12 rows. `n` is 485-1584. Higher in Oct/Nov/Dec

accidents %>% count(hour(timedate))      # 24 rows. `n` is 51-722 (2am, 2pm)

accidents %>% count(gravedad)            # 3 rows. "CON HERIDOS", "CON MUERTOS", "SOLO DANOS"

accidents %>% count(objectid)            # 9982 rows (unique)

accidents %>% count(localidad)           # 19 rows. `n` is 52-1177

accidents %>% count(municipio)           # 1 row ("BOGOTA"). `n` is 9982 --> IGNORE

accidents %>% count(clase)               # 7 rows. `n` is 1 ("INCENDIO") to 8502 ("CHOQUE")

Initial data visialisation, to see if there might be a location/ time pattern

accidents <- accidents %>%
  mutate(week_no = week(timedate), week_day = wday(timedate, label = TRUE))
#1
accidents %>%
  group_by(yday = yday(timedate)) %>%
  summarise(count = n()) %>%
  ggplot(aes(yday, count)) +                           # 'day of the year' (0:1Jan to 365:31Dec)
  geom_point() +
  geom_smooth()

#2
accidents %>%
  filter(month(timedate) == 1) %>%
  group_by(yday = yday(timedate)) %>%
  summarise(count = n()) %>%
  ggplot(aes(yday, count)) +                           # 'day of the year' (0:1Jan to 365:31Dec)
  geom_point() +
  geom_smooth()

#3
ggplot(accidents, aes(timedate)) +
  geom_histogram() +
  facet_wrap(~ localidad)

#4
ggplot(accidents, aes(x, y, colour = localidad)) +     # rough map of the 'localidads'
  geom_point()

Maybe look at a time-based model (especially if there is data for other years too, in which case start with localidad=ENGATIVA). And/or a visualisation (like ACotgreave? if I can get week#). Or time of day (model/viz). Or connected to streetlights (if there’s a joining field)

# Weeks (`week_no`)
accidents %>%
  group_by(week_no, week_day) %>%
  summarise(count = n()) %>%
  ggplot(aes(week_no, week_day)) +
  geom_point(aes(colour = count, size = count, alpha = 0.5))

accidents %>%
  group_by(week_no) %>%
  summarise(count = n()) %>%
  ggplot(aes(week_no, count)) +
  geom_point() +
  geom_smooth()

accidents %>%
  group_by(week_no) %>%
  summarise(count = n()) %>%
  ggplot(aes(week_no, count)) +
  geom_line()

# Hours
#1
accidents %>%
  group_by(week_day, hour = hour(timedate)) %>%
  summarise(count = n()) %>%
  ggplot(aes(hour, week_day)) +
  geom_point(aes(colour = count, size = count, alpha = 0.5))

#2
accidents %>%
  group_by(hour = hour(timedate)) %>%
  summarise(count = n()) %>%
  ggplot(aes(hour)) +
  geom_density()                             

# accidents %>%
#   group_by(hour = hour(timedate)) %>%
#   summarise(count = n()) %>%
#   ggplot(aes(hour)) +
#   geom_density(aes(colour = week_day)) +            # "object 'week_day' not found" ?
#   facet_wrap(~ week_day) # "Error in combine_vars(data, params$plot_env, vars, drop = params$drop) : At least one layer must contain all variables used for facetting" ?"
#3
accidents %>%
  group_by(week_day, hour = hour(timedate)) %>%
  summarise(count = n()) %>%
  ggplot(aes(hour, count)) +
  geom_line(aes(colour = week_day)) +
  facet_wrap(~ week_day)

# Localidads
accidents4 <- accidents %>%
  group_by(month = month(timedate), week_day, hour = hour(timedate), localidad) %>%
  summarise(count = n())
#1
ggplot(accidents4, aes(hour)) +
  geom_density() +
  facet_wrap(~ localidad)

#2
ggplot(accidents4, aes(hour)) +
  geom_density(aes(colour = week_day)) +
  facet_wrap(~ localidad)

#3
ggplot(accidents4, aes(hour)) +
  geom_density() +
  facet_wrap(~ month)

#4
ggplot(accidents4, aes(hour)) +
  geom_density() +
  facet_grid(localidad ~ month)

#5
accidents4 %>%
  filter(localidad == "ENGATIVA") %>%
  ggplot(aes(hour)) +
  geom_density() +
  facet_wrap(~ month)

#6
accidents4 %>%
  filter(localidad == "ENGATIVA") %>%
  ggplot(aes(hour)) +
  geom_density(aes(fill = month)) +
  facet_grid(month ~ .)

Doesn’t appear to be a ‘strong enough’ repeating pattern in 2014, to make it interesting to fit a model to it with modelr. And have been unable to ‘union’ 2014 with the similar-ish 2015&2016 datasets I downloaded, (using Tableau Maestro beta 4), to see if patterns are repeated year-on-year.

So forget about a model and look at the ‘Street Lights’ dataset, to see if that ‘sheds any light’ on (location) patterns in the accidents

lights_raw <- read_csv("data/Street_Lights.csv")

Parsed with column specification:
cols(
  the_geom = col_character(),
  cartodb_id = col_integer(),
  x = col_double(),
  y = col_double(),
  codid = col_integer(),
  direccion = col_character(),
  localidad = col_character(),
  tipointer = col_character(),
  disppdm = col_character(),
  fechaimplementacion = col_date(format = "")
)

lights_raw

#1
ggplot(lights_raw, aes(x, y, colour = localidad)) +     # rough map of the 'localidads'
  geom_point()

#2
ggplot(data=lights_raw, aes(x, y)) +
  geom_point(aes(shape = "."), colour = "red") +
  geom_point(data=accidents, aes(shape = ".", stroke = 0, alpha = 0.5), colour = "black")

#3
accidents_ENGATIVA <- accidents %>%
  filter(localidad == "ENGATIVA")
lights_ENGATIVA <- lights_raw %>%
  filter(localidad == "ENGATIVA")
ggplot(data=lights_ENGATIVA, aes(x, y)) +
  geom_point(aes(shape = "."), size = 3, colour = "black") +
  geom_point(data=accidents_ENGATIVA, aes(shape = ".", stroke = 0, alpha = 0.5), colour = "blue")

#4
deaths_ENGATIVA <- accidents %>%
  filter(localidad == "ENGATIVA") %>%
  filter(gravedad == "CON MUERTOS")
lights_ENGATIVA <- lights_raw %>%
  filter(localidad == "ENGATIVA")
ggplot(data=lights_ENGATIVA, aes(x, y)) +
  geom_point(aes(shape = "."), size = 3, colour = "black") +
  geom_point(data=deaths_ENGATIVA, aes(shape = ".", stroke = 0), size = 2, colour = "red")

There are more accident locations in ENGATIVA than there are street lights. And although some of those street lights have no accidents next to them (including clusters in the SouthEast), the majority does have. So it’s not as simple as saying the patterns are connected (let alone that street lights prevent accidents)

However, only 3 of 11 deaths was next to a street light (in areas where other accidents were happening), so that’s a slightly stronger (visual) case that the lights have a positive effect